Gulf of Maine Sea Surface Temperature Outlook

Region Extent

The spatial extent for Gulf of Maine is displayed below. This bounding box is the same bounding box coordinates used to clip the OISST data when constructing the time series data from the array.

# File paths for various extents based on params$region
poly_path <- switch(
  EXPR = tolower(params$region),
  "gulf of maine"      = "Shapefiles/gmri_sst_focal_areas/apershing_gulf_of_maine.geojson",
  "cpr gulf of maine"  = "Shapefiles/gmri_sst_focal_areas/cpr_gulf_of_maine.geojson",
  "northwest atlantic" = "Shapefiles/gmri_sst_focal_areas/aak_northwest_atlantic.geojson")


# Load the bounding box for Andy's GOM to show they align
poly_path <- paste0(res_path, poly_path)
region_extent <- st_read(poly_path, quiet = TRUE)


# Pull extents for the region for crop
crop_x <- st_bbox(region_extent)[c(1,3)]
crop_y <- st_bbox(region_extent)[c(2,4)]


# Zoom out for cpr extent
if(tolower(params$region) == "cpr gulf of maine"){
  crop_x <- c(-70.875, -65.375)
  crop_y <- c(40.375,   45.125)}


# Full plot
ggplot() +
  geom_sf(data = new_england, fill = "gray90") +
  geom_sf(data = canada, fill = "gray90") +
  geom_sf(data = greenland, fill = "gray90") +
  geom_sf(data = region_extent, color = "royalblue", 
          fill = "transparent", linetype = 2, size = 1) +
  map_theme +
  coord_sf(xlim = crop_x, 
           ylim = crop_y, expand = T)

GMRI OISST Resources on Box

For many of our most-used areas, regional climatologies and anomalies have been processed using a handful of jupyter notebook work flows to take advantage of {xarray}. The following sections will highlight what those work flows are, what they detail, and how they can be updated using the Gulf of Maine as an example.

The root folder for the following OISST products is ~/Box/RES_Data/OISST/oisst_mainstays.

Regional Timeseries

A region-masked time series is the most basic building block for relaying information for any of these areas. For any desired area (represented by a spatial polygon) a time series table of observed sea surface temperature is returned that contains the mean sea surface temperature for that area only for each time step.

Rather than repeat this process again for the climatology table and the anomalies these three values and a few additional fields are stored in the same table.

Pre-processed tables have been placed on box as part of the NSF C-ACCEL project under the folder: ~/Box/RES_Data/OISST/oisst_mainstays/regional_timeseries. A function to quickly access these timelines has been added to the {gmRi} package as oisst_access_timeseries().

# Switch path based on local vs. docker
timeseries_path <- switch(
  EXPR = tolower(params$region),
  "gulf of maine"      = "OISSTv2_anom_apershing_gulf_of_maine.csv",
  "cpr gulf of maine"  = "OISSTv2_anom_cpr_gulf_of_maine.csv",
  "northwest atlantic" = "OISSTv2_anom_northwest_atlantic.csv")

# Append to form complete path
timeseries_path <- str_c(
  oisst_path, 
  "regional_timeseries/gmri_sst_focal_areas/", 
  timeseries_path)
  

# Load timeseries
region_timeseries <- read_csv(timeseries_path,
                              guess_max = 1e5,
                              col_types = cols())

# Display Table of first 6 entries
head(region_timeseries) %>% 
  mutate_if(is_numeric, round, 2) %>% 
  kable() %>% 
  kable_styling()
time modified_ordinal_day sst sst_clim sst_anom clim_sd log_lik
1981-09-01 245 15.78 16.72 -0.94 2.32 1.84
1981-09-02 246 15.79 16.67 -0.89 2.29 1.82
1981-09-03 247 15.49 16.64 -1.14 2.28 1.87
1981-09-04 248 14.99 16.60 -1.60 2.25 1.98
1981-09-05 249 14.84 16.51 -1.67 2.21 2.00
1981-09-06 250 15.08 16.49 -1.41 2.18 1.91

Climatologies are currently set up to calculate daily averages on a modified julian day, such that every March 1st and all days after fall on the same day, regardless of whether it is a leap year or not. This preserves comparisons across calendar dates such-as: “The average temperature on march 1st is ____ for the reference period ____ to ____”

In these tables sst is the mean temperature observed for that date averaged across all cells within the area. sst_clim & clim_sd are the climate means and standard deviations for a 1982-2011 climatology. sst_anom is the daily observed minus the climate mean, and log_lik is the log likelihood of getting that value given a normal distribution with mean of sst_clim and standard deviation of clim_sd.

Warming Rates

Regional warming trends below were calculated using all the available data from Fall of 1981 through the end of 2020.

Annual

# Summarise by year to return mean annual anomalies and variance
annual_summary <- region_timeseries %>% 
  mutate(year = year(time)) %>% 
  group_by(year) %>% 
  summarise(sst = mean(sst, na.rm = T),
            sst_anom = mean(sst_anom, na.rm = T), 
            .groups = "keep") %>% ungroup() %>% 
  mutate(yr_as_dtime = as.Date(paste0(year, "-07-02")))


# Build Regression Equations
lm_all <- lm(sst_anom ~ year, data = annual_summary) %>% 
  coef() %>% 
  round(3)
lm_15  <- lm(sst_anom ~ year, data = filter(annual_summary, year %in% c(2006:2020))) %>% 
  coef() %>% 
  round(3)

# Equation to paste in
eq_all <- paste0("y = ", lm_all['(Intercept)'], " + ", lm_all['year'], " x")
eq_15  <- paste0("y = ", lm_15['(Intercept)'], " + ", lm_15['year'], " x")


# Plot the warming rate for the region
ggplot(data = annual_summary, aes(yr_as_dtime, sst_anom)) +
  
  # Add daily data
  geom_line(data = region_timeseries,
            aes(time, sst_anom),
            alpha = 0.5, color = "gray") +
  
  # Overlay yearly means
  geom_line(alpha = 0.7) +
  geom_point(alpha = 0.7) +
  
  # Add regression lines 
  geom_smooth(method = "lm",
              aes(color = "1982-2020 Regional Trend"),
              formula = y ~ x, se = F,
              linetype = 2) +
  geom_smooth(data = filter(annual_summary, year %in% c(2006:2020)),
              method = "lm", 
              aes(color = "2006-2020 Regional Trend"),
              formula = y ~ x, se = F, 
              linetype = 2) +
  # Manually add equations so they show yearly not daily coeff
  geom_text(data=data.frame(), aes(label = eq_all, x = min(region_timeseries$time), y = Inf),
            hjust = 0, vjust = 2, color = gmri_cols("gmri blue")) +
  geom_text(data=data.frame(), aes(label = eq_15, x = min(region_timeseries$time), y = Inf),
            hjust = 0, vjust = 3.5, color = gmri_cols("orange")) +

  # Colors, labels, theme
  scale_color_manual(values = c(
    "1982-2020 Regional Trend" = as.character(gmri_cols("gmri blue")),
    "2006-2020 Regional Trend" = as.character(gmri_cols("orange")))) +
  labs(x = "", y = "Sea Surface Temperature Anomaly",
       caption = expression("Regression Coefficients Reflect Annual Change in"~degree~C)) +
  theme(legend.title = element_blank(),
        legend.position = c(0.85, 0.1),
        legend.background = element_rect(fill = "transparent"),
        legend.key = element_rect(fill = "transparent", color = "transparent"))

# Removed when daily data was added
# # Add equation for full timeseries using the input data
# stat_poly_eq(formula = y ~ x,
#              aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~~")),
#              color = gmri_cols("gmri blue"),
#              parse = T) +
# Add regression line - all data
# Add equation for last 15 years
# stat_poly_eq(formula = y ~ x,
#              aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~~~")),
#              color = gmri_cols("orange"),
#              parse = T, label.x = 0.05, label.y = 0.9) +

Quarterly

# Repeat for the seasons (equal quarters here*)
quarter_summary <- region_timeseries %>% 
  mutate(year = year(time),
         season = factor(quarter(time, fiscal_start = 1)),
         season = fct_recode(season, c("Jan 1 - March 31" = "1"), c("Apr 1 - Jun 30" = "2"), 
                             c("Jul 1 - Sep 30" = "3"), c("Oct 1 - Dec 31" = "4"))) %>% 
  group_by(year, season) %>% 
  summarise(sst = mean(sst, na.rm = T),
            sst_anom = mean(sst_anom, na.rm = T), 
            .groups = "keep") 

# Plot
quarter_summary %>% 
  ggplot(aes(year, sst_anom)) +
  geom_line(group = 1) +
  geom_point() +
  geom_smooth(method = "lm", 
              aes(color = "Regional Trend"),
              formula = y ~ x, se = F, linetype = 2) +
  stat_poly_eq(formula = y ~ x,
               color = gmri_cols("gmri blue"),
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
               parse = T) +
  scale_color_manual(values = c("Regional Trend" = as.character(gmri_cols("orange")))) +
  labs(x = "", y = "Sea Surface Temperature Anomaly") +
  theme(legend.title = element_blank(),
        legend.position = c(0.875, 0.05),
        legend.background = element_rect(fill = "transparent"),
        legend.key = element_rect(fill = "transparent", color = "transparent")) +
  facet_wrap(~season)

Marine Heatwaves

The {heatwaveR} package provides a relatively quick way of working with tabular data to calculate a seasonal climate mean, as well as heatwaves and coldwaves at a desired threshold.

These functions can be used to get the climatology using the standard day of year if desired, and track the number and length of heatwave events. Heatwave events follow the definition of Hobday et al. 2016 and are set up to use a 90% threshold for heatwave/coldwave detection.

The heatwave history for the Gulf of Maine trawl region displayed above is as follows:

# Wrapper function to do heatwaves and coldwaves simultaneously at 90%
pull_heatwave_events <- function(temperature_timeseries) {
    
    # Pull the two column dataframe for mhw estimation
    test_ts <- data.frame(t = temperature_timeseries$time, 
                          temp = temperature_timeseries$sst)
    
    # Detect the events in a time series
    ts  <- ts2clm(data = test_ts, climatologyPeriod = c("1982-01-01", "2011-12-31"))
    
    #heatwaves
    mhw <- detect_event(ts)                         
    
    # prep heatwave data
    mhw_out <- mhw$climatology %>% 
        mutate(sst_anom = temp - seas) %>% 
        select(time = t,
               sst = temp,
               seas,
               sst_anom,
               mhw_thresh = thresh,
               mhw_event = event,
               mhw_event_no = event_no)

    # 2. Detect cold spells
    ts <- ts2clm(data = test_ts, 
                 climatologyPeriod = c("1982-01-01", "2011-12-31"), 
                 pctile = 10)
    mcs <- detect_event(ts, coldSpells = TRUE) #coldSpells = TRUE flips boolean to < thresh
    
    # prep cold spell data
    mcs_out <- mcs$climatology %>%
        select(time = t,
               mcs_thresh = thresh,
               mcs_event = event,
               mcs_event_no = event_no)
    

    # join heatwaves to coldwaves
    hot_and_cold <- left_join(mhw_out, mcs_out, by = "time")
    
    
    
    # 3. Data formatting for plotting, 
    # adds columns to plot hw and cs seperately
    events_out <- hot_and_cold %>% 
        mutate(status = ifelse(mhw_event == TRUE, "Marine Heatwave Event", "Sea Surface Temperature"),
               status = ifelse(mcs_event == TRUE, "Marine Cold Spell Event", status),
               hwe = ifelse(mhw_event == TRUE, sst, NA),
               cse = ifelse(mcs_event == TRUE, sst, NA),
               nonevent = ifelse(mhw_event == FALSE & mcs_event == FALSE, sst, NA)) 
    
    # Close the gaps between a mhw event and sst (might not need if full line for temp exists)
    events_out <- events_out %>% 
        mutate(hwe = ifelse(is.na(hwe) & is.na(lag(hwe)) == FALSE, sst, hwe),
               cse = ifelse(is.na(cse) & is.na(lag(cse)) == FALSE, sst, cse))
    
    
    return(events_out)
}


# Use function to process heatwave data for plotting
region_heatwaves <- pull_heatwave_events(region_timeseries)

Heatwave timelines can be then be plotted using the {ggplot2} package for static plots, or by using the {plotly} package for interactive content.

Interactive SST Timeline

For anything we wish to host on the web there is an option to display tables and graphs that are interactive. The {plotly} package is one-such tool for producing plots that allow users to pan, zoom, and highlight discrete observations.

# Helper function
plotly_mhw_plots <- function(data){
    
    # How to get rgb() from colors
    plot_cols <- list(
        "gray20"    = col2rgb(col = "gray20"),
        "gray40"    = col2rgb(col = "gray40"),
        "royalblue" = col2rgb(col = "royalblue"),
        "darkred"   = col2rgb(col = "darkred"),
        "royalblue" = col2rgb(col = "royalblue"))
    
    
    # Building the plot
    fig <- plot_ly(data, x = ~time) 
    
    # Sea Surface Temperature
    fig <- fig %>% add_trace(y = ~sst, 
                             name = 'Sea Surface Temperature',
                             mode = 'lines', 
                             type = "scatter",
                             line = list(color = "rgb(65, 105, 225)", 
                                         width = 2)) 
    # Heatwave Threshold
    fig <- fig %>% add_trace(y = ~mhw_thresh, 
                             name = 'MHW Threshold ', 
                             mode = 'lines', 
                             type = "scatter",
                             line = list(color = "rgb(51, 51, 51)", 
                                         width = 1, 
                                         dash = 'dot')) 
    # Seasonal Climatology
    fig <- fig %>% add_trace(y = ~seas, 
                             name = 'Seasonal Climatology ', 
                             mode = 'lines', 
                             type = "scatter",
                             line = list(color = "rgb(102, 102, 102)", 
                                         width = 3, 
                                         dash = 'dash')) 
    # Marine Cold Spell Threshold
    fig <- fig %>% add_trace(y = ~mcs_thresh, 
                             name = 'MCS Threshold ', 
                             mode = 'lines', 
                             type = "scatter",
                             line = list(color = "rgb(51, 51, 51)", 
                                         width = 1, 
                                         dash = 'dot')) 
    # Heatwave Event
    fig <- fig %>% add_trace(y = ~hwe, 
                             name = 'Marine Heatwave Event ', 
                             mode = 'lines', 
                             type = "scatter",
                             line = list(color = "rgb(139, 0, 0)", 
                                         width = 2)) 
   
    # Cold Spell Event
    fig <- fig %>% add_trace(y = ~cse, 
                             name = 'Marine Cold Spell Event ', 
                             mode = 'lines', 
                             type = "scatter",
                             line = list(color = "rgb(65, 105, 225)", 
                                         width = 2)) 
    
    # Axis Formatting
    fig <- fig %>% layout(xaxis = list(title = ""),
                          yaxis = list (title = "Temperature (degrees C)"))
    
    
    # Legend formatting
    fig <- fig %>% layout(legend = list(orientation = 'h'))
    
    
    return(fig)
}



# Grab data from the most recent year through present day
last_year <- Sys.Date() - 365
last_year <- last_year - yday(last_year) + 1
last_yr_heatwaves <- region_heatwaves %>% 
  filter(time >= last_year)

# Plot that timeseries
last_yr_heatwaves %>% 
  filter(time >= last_year) %>% 
  plotly_mhw_plots()

Static SST Timeline

In most cases a static image will suffice, in which case {ggplot2} is a familiar go-to.

# Set colors by name
color_vals <- c(
  "Sea Surface Temperature" = "royalblue",
  "Heatwave Event"          = "darkred",
  "Cold Spell Event"        = "lightblue",
  "MHW Threshold"           = "gray30",
  "MCS Threshold"           = "gray30",
  "Seasonal Climatology"    = "gray30")

# Set the label with degree symbol
ylab <- expression("Sea Surface Temperature"~degree~C)


# How many heatwave events:
num_events <- max(last_yr_heatwaves$mhw_event_no, na.rm = T) - min(last_yr_heatwaves$mhw_event_no, na.rm = T)

# How many heatwave days
num_hw_days <- sum(last_yr_heatwaves$mhw_event, na.rm = T)


# Plot the last 365 days
ggplot(last_yr_heatwaves, aes(x = time)) +
    geom_segment(aes(x = time, xend = time, y = seas, yend = sst), 
                 color = "royalblue", alpha = 0.25) +
    geom_segment(aes(x = time, xend = time, y = mhw_thresh, yend = hwe), 
                 color = "darkred", alpha = 0.25) +
    geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
    geom_line(aes(y = hwe, color = "Heatwave Event")) +
    geom_line(aes(y = cse, color = "Cold Spell Event")) +
    geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 2, size = .5) +
    geom_line(aes(y = mcs_thresh, color = "MCS Threshold"), lty = 2, size = .5) +
    geom_line(aes(y = seas, color = "Seasonal Climatology"), size = .5) +
    scale_color_manual(values = color_vals) +
    scale_x_date(date_labels = "%b", date_breaks = "1 month") +
    theme(legend.title = element_blank(),
          legend.position = "bottom") +
    labs(x = "", y = ylab, 
         caption = paste0("X-axis start date: ", last_year, 
                          "\nNumber of Heatwave events: ", num_events,
                          "\nNumber of Heatwave Days: ", num_hw_days))

Heatwave Events

As mentioned above heatwave events were determined using the methods of Hobday et al. 2016 and implemented using {heatwaveR} with a 90% threshold. The region used is again that same polygon outlined in the first figure, the bottom right panel of the maps at the beginning of the document.

# Prep the legend title
guide_lab <- expression("Sea Surface Temperature Anomaly"~degree~C)

# get new axis dimensions, y = year, x = day within year
# use flate_date so that they don't stair step
base_date <- as.Date("2000-01-01")
grid_data <- region_heatwaves %>% 
  mutate(year = year(time),
         yday = yday(time),
         flat_date = as.Date(yday-1, origin = base_date))


# assemble plot
ggplot(grid_data, aes(x = flat_date, y = year)) +
  
  # background box fill for missing dates
  geom_rect(xmin = base_date, xmax = base_date + 365, 
            ymin = min(grid_data$year) - .5, ymax = max(grid_data$year) + .5, 
            fill = "gray75", color = "transparent") +
  
  # tile for sst colors
  geom_tile(aes(fill = sst_anom)) +
  
  # points for heatwave events
  geom_point(data = filter(grid_data, mhw_event == TRUE),
             aes(x = flat_date, y = year), size = .25)  +
  scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
  scale_y_continuous(limits = c(1980.5, 2020.5), expand = c(0,0)) +
  labs(x = "", y = "") +
  
  #scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
  scale_fill_distiller(palette = "RdBu", na.value = "transparent") +
  
  #5 inches is default rmarkdown height for barheight
  guides("fill" = guide_colorbar(title = guide_lab, 
                                 title.position = "right", 
                                 title.hjust = 0.5,
                                 barheight = unit(4.8, "inches"), 
                                 frame.colour = "black", 
                                 ticks.colour = "black")) +  
  theme_classic() +
  theme(legend.title = element_text(angle = 90))

Anomaly Maps

For demonstration purposes I will only be loading the 2020 global anomaly data for the below images.

# Access information to netcdf on box
nc_year <- "2020"
anom_path <- str_c(oisst_path, "annual_anomalies/1982to2011_climatology/daily_anoms_", nc_year, ".nc")


# Load 2020 as stack
anoms_2020 <- stack(anom_path)

Mean Global

The following code will subset the anomalies for July and plot the average sea surface temperature anomalies for that month:

# Get the mean temperature anomalies for July
july_dates <- which(str_sub(names(anoms_2020), 7, 8) == "07")
july_avg <- mean(anoms_2020[[july_dates]])

# Convert wgs84 raster to stars array
july_st <- st_as_stars(rotate(july_avg))



# Plot global map
ggplot() +
  geom_stars(data = july_st) +
  geom_sf(data = world, fill = "gray90") +
  # scale_fill_gradient2(low = "blue",
  #                      mid = "white",
  #                      high = "red") +
  scale_fill_distiller(palette = "RdBu", na.value = "transparent") +
  map_theme +
  coord_sf(expand = FALSE) +
  guides("fill" = guide_colorbar(title = "Average Sea Surface Temperature Anomaly",
                                 title.position = "top", 
                                 title.hjust = 0.5,
                                 barwidth = unit(4, "in"), 
                                 frame.colour = "black", 
                                 ticks.colour = "black"))

Regional Map

Using the July average sst anomaly we can then clip down to the Gulf of Maine for an aerial view of the region. For this execution the bounding box of lat-lon coordinates from the regional extent will be used.

# Clip Raster - Convert to stars
shape_extent <- c(crop_x, crop_y)
region_ras <- crop(rotate(july_avg), extent(shape_extent))
region_st <- st_as_stars(region_ras)

# Get crop bounds for coord_sf
crop_x <- st_bbox(region_st)[c(1,3)] 
crop_y <- st_bbox(region_st)[c(2,4)]


# Zoom out for cpr extent, same as Andy's GOM
if(tolower(params$region) == "cpr gulf of maine"){
  crop_x <- c(-71, -65.5)
  crop_y <- c(40.5, 45)}


# Plot Gulf of Maine - wgs84
ggplot() +
  geom_stars(data = region_st) +
  geom_sf(data = new_england, fill = "gray90") +
  geom_sf(data = canada, fill = "gray90") +
  geom_sf(data = greenland, fill = "gray90") +
  scale_fill_distiller(palette = "RdBu", na.value = "transparent") +
  map_theme +
  coord_sf(xlim = crop_x, ylim = crop_y, expand = T) +
  guides("fill" = guide_colorbar(
    title = "Average Sea Surface Temperature Anomaly",
    title.position = "top",
    title.hjust = 0.5,
    barwidth = unit(4, "in"),
    frame.colour = "black",
    ticks.colour = "black"))

Regionally Projected Map

These lat/lon region extent boxes are good showcases of how something that is “rectangular” by its lat/lon dimensions that has its shape change to when projected onto a non-flat surface of the earth.

To get the OISST anomaly grid to display on a curviliniear projection the originals were cropped using the corner coordinates of the extent plotted above. The raster was then warped into a regular grid using an appropriate projection as displayed below.

For the most dramatic case set the region parameter to “Northwest Atlantic”

#### Candidate coordinate reference systems ####

####  NOTE: 
# transformation from rasters different than polygons 
# as cells need to stretch and bend  ####
# Resource: Tranformation vs. Warping 
# https://r-spatial.github.io/stars/articles/stars5.html

# Robinson projection
robinson_proj <- "+proj=robin"

# Albers equal area: centered on -70 degrees
# The settings for lat_1 and lat_2 are the locations at which 
# the cone intersects the earth, so distortion is minimized at those latitudes
alb_70 <- "+proj=aea +lat_1=30 +lat_2=50 +lon_0=-70"

# custom sterographic projection for this nw atlantic, centered using lon_0
stereographic_north <- "+proj=stere +lat_0=90 +lat_ts=75 +lon_0=-57"

# equal earth projection - not working, don't think sf has functionality for it yet
eqearth_proj <- "+proj=eqearth"

# Choose crs using parameter
projection_crs <- switch(
  tolower(params$region),
  "gulf of maine" = alb_70,
  "cpr gulf of maine" = alb_70,
  "northwest atlantic" = stereographic_north)


# Transform all the polygons
region_projected        <- st_transform(region_extent, crs = projection_crs)
canada_projected        <- st_transform(canada, crs = projection_crs)
newengland_projected    <- st_transform(new_england, crs = projection_crs)
greenland_projected     <- st_transform(greenland, crs = projection_crs)

# coord_sf Crop bounds in projection units for coord_sf
crop_x <- st_bbox(region_projected)[c(1,3)] 
crop_y <- st_bbox(region_projected)[c(2,4)]


# Zoom out for cpr extent, same as Andy's GOM
if(tolower(params$region) == "cpr gulf of maine"){
  crop_x <- c(-73185.13, 386673.34)
  crop_y <- c(4269044,   4813146)}

# Lower the ymin a touch for NW Atlantic
if(tolower(params$region) == "northwest atlantic"){crop_y <- crop_y - c(100000, 0)}


# Warp to grid of chosen CRS
projection_grid <- region_st %>% 
  st_transform(projection_crs) %>% 
  st_bbox() %>%
  st_as_stars()
region_warp_ras <- region_st %>% 
  st_warp(projection_grid) 


# Plot everything together
ggplot() +
  geom_stars(data = region_warp_ras) +
  geom_sf(data = newengland_projected, fill = "gray90") +
  geom_sf(data = canada_projected, fill = "gray90") +
  geom_sf(data = greenland_projected, fill = "gray90") +
  coord_sf(crs = projection_crs, 
           xlim = crop_x, ylim = crop_y, expand = T) +
  map_theme +
  scale_fill_distiller(palette = "RdBu", na.value = "transparent") +
  guides("fill" = guide_colorbar(title = "Average Sea Surface Temperature Anomaly",
                                 title.position = "top", 
                                 title.hjust = 0.5,
                                 barwidth = unit(4, "in"), 
                                 frame.colour = "black", 
                                 ticks.colour = "black")) +
  theme(
    panel.border = element_rect(color = "black", fill = NA),
    plot.background = element_rect(color = "transparent", fill = "transparent"),
    axis.title.x = element_blank(), # turn off titles
    axis.title.y = element_blank(),
    legend.position = "bottom", 
    legend.title.align = 0.5)

Warming Rates

Global warming rates are calculated by getting the average temperature for each year across all cells in the data, then using a linear regression for each to determine the annual warming rate in degrees Celsius.

The warming rates are then ranked from high to low to identify areas with the most rapid ocean warming. The following maps display the warming rates and the percentile ranks for the top 20% of warming rates, or rates above the 80th percentile of the data.

# Access information to warming rate netcdf files on box
rates_path <- paste0(oisst_path, "warming_rates/annual_warming_rates")


# Going to reclassify raster from continuous to discrete: 
# Source: https://www.earthdatascience.org/courses/earth-analytics/lidar-raster-data-r/classify-raster/

reclassify_to_discrete <- function(ranks_stack, rates_stack, percentile_cutoff = 80){

  # create classification matrix
  reclass_df <- c(0.00, 0.05,  0, #####
                  0.05, 0.10,  5,
                  0.10, 0.15, 10,
                  0.15, 0.20, 15,
                  0.20, 0.25, 20,
                  0.25, 0.30, 25,
                  0.30, 0.35, 30,
                  0.35, 0.40, 35,
                  0.40, 0.45, 40,
                  0.45, 0.50, 45,
                  0.50, 0.55, 50,
                  0.55, 0.60, 55,
                  0.60, 0.65, 60,
                  0.65, 0.70, 65,
                  0.70, 0.75, 70,
                  0.75, 0.80, 75,
                  0.80, 0.85, 80,
                  0.85, 0.90, 85,
                  0.90, 0.95, 90,
                  0.95,    1, 95)
  #####
  
  
  # reshape the object into a matrix with columns and rows
  reclass_m <- matrix(reclass_df,
                  ncol = 3,
                  byrow = TRUE)
  
  # reclassify the raster using the reclass object - reclass_m
  ranks_classified <- reclassify(ranks_stack,
                                 reclass_m)
  
  # Save the un-masked layers as stars objects
  rates_raw_st <- st_as_stars(rotate(rates_stack))
  ranks_raw_st <- st_as_stars(rotate(ranks_stack))
  
  # Masking Below percentile cutoff - for ranking raster and rates raster
  ranks_stack[ranks_classified < percentile_cutoff] <- NA
  rates_stack[ranks_classified < percentile_cutoff] <- NA
  
  # Converting to stars
  rates_st   <- st_as_stars(rotate(rates_stack))
  ranks_st   <- st_as_stars(rotate(ranks_stack))
  ranks_c_st <- st_as_stars(rotate(ranks_classified))
  
  # #get scale ranges so they still pop
  # rate_range <- c("min" = cellStats(rates_stack, 'min')[1], 
  #                 "max" = cellStats(rates_stack,'max')[1])
  
  # Return the three stars objects
  return(list("rates_raw" = rates_raw_st,
              "ranks_raw" = ranks_raw_st,
              "rates" = rates_st, 
              "ranks" = ranks_st, 
              "ranks_discrete" = ranks_c_st))

}

1982-2011

####  Data Prep  ####
percentile_cutoff <- 80

# 1982-2020
rates_stack_all <- stack(str_c(rates_path, "1982to2020.nc"), 
                         varname = "annual_warming_rate")
ranks_stack_all <- stack(str_c(rates_path, "1982to2020.nc"), 
                         varname = "rate_percentile")

rates_prepped <- reclassify_to_discrete(rates_stack = rates_stack_all, 
                                        ranks_stack = ranks_stack_all,
                                        percentile_cutoff = percentile_cutoff)

rates_st <- rates_prepped$rates
rates_raw_st <- rates_prepped$rates_raw
ranks_st <- rates_prepped$ranks

Global Rates

####  Plot Setup  ####

# guide label
guide_lab <- expression("Annual Warming Rate - "~degree~C)

# ggplot using stars
ggplot() +
  geom_stars(data = rates_st) +
  #geom_stars(data = rates_raw_st) +
  geom_sf(data = world, fill = "gray90") + 
  scale_fill_viridis_c(na.value = "black", option = "plasma") +
  guides("fill" = guide_colorbar(title = guide_lab,
                                 title.position = "top", 
                                 title.hjust = 0.5,
                                 barwidth = unit(4, "in"), 
                                 frame.colour = "black", 
                                 ticks.colour = "black")) +
  map_theme +
  coord_sf(expand = FALSE) +
  labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown.",
                       "\n Data from Years: 1982-2020"))

Global Ranks

# ggplot using stars
ggplot() +
  geom_stars(data = ranks_st) +
  geom_sf(data = world, fill = "gray90") + 
  scale_fill_viridis_c(na.value = "black", option = "plasma") +
  guides("fill" = guide_colorbar(title = "Percentile of Annual Warming Rate",
                                 title.position = "top", 
                                 title.hjust = 0.5,
                                 barwidth = unit(4, "in"), 
                                 frame.colour = "black", 
                                 ticks.colour = "black")) +
  map_theme +
  coord_sf(expand = FALSE) +
  labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown.",
                       "\n Data from Years: 1982-2020"))

2006-2020

####  Data Prep  ####

# 2006-2020
rates_stack_all <- stack(str_c(rates_path, "2006to2020.nc"), 
                         varname = "annual_warming_rate")
ranks_stack_all <- stack(str_c(rates_path, "2006to2020.nc"), 
                         varname = "rate_percentile")

rates_prepped <- reclassify_to_discrete(rates_stack = rates_stack_all, 
                                        ranks_stack = ranks_stack_all)

rates_st <- rates_prepped$rates
ranks_st <- rates_prepped$ranks

Global Rates

####  Plot Setup  ####

# guide label
guide_lab <- expression("Annual Warming Rate - "~degree~C)

# ggplot using stars
ggplot() +
  geom_stars(data = rates_st) +
  geom_sf(data = world, fill = "gray90") + 
  scale_fill_viridis_c(na.value = "black", option = "plasma") +
  guides("fill" = guide_colorbar(title = guide_lab,
                                 title.position = "top", 
                                 title.hjust = 0.5,
                                 barwidth = unit(4, "in"), 
                                 frame.colour = "black", 
                                 ticks.colour = "black")) +
  map_theme +
  coord_sf(expand = FALSE) +
  labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown.",
                       "\n Data from Years: 2006-2020"))

Global Ranks

# ggplot using stars
ggplot() +
  geom_stars(data = ranks_st) +
  geom_sf(data = world, fill = "gray90") + 
  scale_fill_viridis_c(na.value = "black", option = "plasma") +
  guides("fill" = guide_colorbar(title = "Percentile of Annual Warming Rate",
                                 title.position = "top", 
                                 title.hjust = 0.5,
                                 barwidth = unit(4, "in"), 
                                 frame.colour = "black", 
                                 ticks.colour = "black")) +
  map_theme +
  coord_sf(expand = FALSE) +
  labs(caption = str_c("Top ", percentile_cutoff, "% of Warming Rates Shown.",
                       "\n Data from Years: 2006-2020"))

 

A work by Adam A. Kemberling

Akemberling@gmri.org